//***********************************************************//
//****      COMPUTE CONTRIBUTION OF NEW ENTREPRENEURS     ***//
//***********************************************************//

/*
    This file computes the contribution to different indicator of newly entrepreneurs.
    We compute survival rate, contribution to production, capital accumulation etc.
*/

// RECALL //
// GROUP1 = counterfactual group -> same fraction as in the benchmark model //
// GROUP2 = include selection of new entrepreneurs //
// GROUP3 = total fraction //




// INPUT:: distEgroup1, distEgroup2, distEgroup3, distEgroup4
#if NOPOLICY == 1
void SIMULATION_ENT(double *newJNE1, double *newJE1, double *newJNE2, double *newJE2, double *newJNE3, double *newJE3, double *valueJNE, double *valueJE, double *valueUiE, double *valueUnE, double *valueWWE, double *valueWWNE, double *valueUnNE, double *valueUiNE, double *defaultJNE, double *collateralJNE, double *sloanJNE, double *saveresJNE, double *searchJNE_W, int *isaveJNE, double *collateralJE, double *saveresJE, double *searchJE_W, int *isaveJE, double aggcapital, double agglabor, int type_simul) {
#endif

#if SEA == 1
void SIMULATION_ENT(double *newJNE1, double *newJE1, double *newJNE2, double *newJE2, double *newJNE3, double *newJE3, double *valueJNE, double *valueJE, double *valueUiE, double *valueUnE, double *valueWWE, double *valueWWNE, double *valueUnNE, double *valueUiNE,  double *defaultJNE, double *collateralJNE, double *sloanJNE, double *saveresJNE, double *searchJNE_W, int *isaveJNE, double *collateralJE, double *saveresJE, double *searchJE_W, int *isaveJE, double *defaultJNEi, double *collateralJNEi, double *sloanJNEi, double *saveresJNEi, double *searchJNEi_W, int *isaveJNEi, double *collateralJEi, double *saveresJEi, double *searchJEi_W, int *isaveJEi, double *valueJNEi, double *valueJEi, double *rateJNEi, double aggcapital, double agglabor, int type_simul) {
#endif




printf("Initialization...\n");

double *zeroE, expprob;  // vector of zeros (at each new ii, ee, tt), reinitialize
zeroE = (double *) calloc((ifulldimETsimul), sizeof(double)); // calloc already allocates zero values


// READINPUT DISTEGROUP1,2,3,4...
// INSURED GROUPS //
// GROUP 1 = entrepreneur would haven't entered entrepreneurship without insurance (purely selected entrepreneurs)
double *distEgroup1JNEi, *distEgroup1JEi, *distEgroup1JNE, *distEgroup1JE;
distEgroup1JNEi = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup1JEi = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup1JNE = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup1JE = (double *) calloc((ifulldimETsimul), sizeof(double));

// GROUP 2 = entrepreneur would have entered entrepreneurship even without insurance (control group)
double *distEgroup2JNEi, *distEgroup2JEi, *distEgroup2JNE, *distEgroup2JE;
distEgroup2JNEi = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup2JEi = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup2JNE = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup2JE = (double *) calloc((ifulldimETsimul), sizeof(double));

// UNINSURED GROUPS //
// GROUP 3 = entrepreneurs entering entrepreneurship in total
double *distEgroup3JNE, *distEgroup3JE;
distEgroup3JNE = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup3JE = (double *) calloc((ifulldimETsimul), sizeof(double));

// GROUP 4 = entrepreneur entering entrepreneurship as Ui.
double *distEgroup4JNE, *distEgroup4JE;
distEgroup4JNE = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup4JE = (double *) calloc((ifulldimETsimul), sizeof(double));

// TEMP VALUES //
// INSURED GROUPS //
// GROUP 1 = entrepreneur would haven't entered entrepreneurship without insurance (purely selected entrepreneurs)
double *distEgroup1JNEitemp, *distEgroup1JEitemp, *distEgroup1JNEtemp, *distEgroup1JEtemp;
distEgroup1JNEitemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup1JEitemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup1JNEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup1JEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));

// GROUP 2 = entrepreneur would have entered entrepreneurship even without insurance (control group)
double *distEgroup2JNEitemp, *distEgroup2JEitemp, *distEgroup2JNEtemp, *distEgroup2JEtemp;
distEgroup2JNEitemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup2JEitemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup2JNEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup2JEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));

// UNINSURED GROUPS //
// GROUP 3 = entrepreneurs entering entrepreneurship in total
double *distEgroup3JNEtemp, *distEgroup3JEtemp;
distEgroup3JNEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup3JEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));

// GROUP 4 = entrepreneur entering entrepreneurship as Ui.
double *distEgroup4JNEtemp, *distEgroup4JEtemp;
distEgroup4JNEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));
distEgroup4JEtemp = (double *) calloc((ifulldimETsimul), sizeof(double));



// INITIALIZE STATISTICS //
    double avgW_group1, avgW_group2, avgW_group3, avgW_group4, distgroup1, distgroup2, distgroup3, distgroup4, avgT_group1, avgT_group2, avgT_group3, avgT_group4, avgY_group1, avgY_group2, avgY_group3, avgY_group4, avgK_group1, avgK_group2, avgK_group3, avgK_group4, avgMPLDL_group1, avgMPLDL_group2, avgMPLDL_group3, avgMPLDL_group4, avgY_group1_bef, avgY_group2_bef, avgY_group3_bef, avgY_group4_bef, avgK_group1_bef, avgK_group2_bef, avgK_group3_bef, avgK_group4_bef, avgBANKRUPT_group1, avgBANKRUPT_group2, avgBANKRUPT_group3, avgBANKRUPT_group4;

double dist1[maxtimesimul], dist2[maxtimesimul], dist3[maxtimesimul], dist4[maxtimesimul];
double dist1T[maxtimesimul*maxtheta], dist2T[maxtimesimul*maxtheta], dist3T[maxtimesimul*maxtheta], dist4T[maxtimesimul*maxtheta];


double *avgW_g1, *avgW_g2, *avgW_g3, *avgW_g4, *dist_group1, *dist_group2, *dist_group3, *dist_group4, *avgT_g1, *avgT_g2, *avgT_g3, *avgT_g4, *avgY_g1, *avgY_g2, *avgY_g3, *avgY_g4, *avgBANKRUPT_g1, *avgBANKRUPT_g2, *avgBANKRUPT_g3, *avgBANKRUPT_g4, *avgK_g1, *avgK_g2, *avgK_g3, *avgK_g4, *avgMPLDL_g1, *avgMPLDL_g2, *avgMPLDL_g3, *avgMPLDL_g4, *growthK_group1, *growthK_group2, *growthK_group3, *growthK_group4, *growthY_group1, *growthY_group2, *growthY_group3, *growthY_group4, *surv_group1, *surv_group2, *surv_group3, *surv_group4;

// avg wealth level //
avgW_g1 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgW_g2 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgW_g3 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgW_g4 = (double *) calloc((ifulldimETsimul), sizeof(double));

// Total dist //
dist_group1 = (double *) calloc((ifulldimETsimul), sizeof(double)); dist_group2 = (double *) calloc((ifulldimETsimul), sizeof(double)); dist_group3 = (double *) calloc((ifulldimETsimul), sizeof(double)); dist_group4 = (double *) calloc((ifulldimETsimul), sizeof(double));

// Average skill //
avgT_g1 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgT_g2 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgT_g3 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgT_g4 = (double *) calloc((ifulldimETsimul), sizeof(double));
    
// Average prod // OK
avgY_g1 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgY_g2 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgY_g3 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgY_g4 = (double *) calloc((ifulldimETsimul), sizeof(double));
    
// avgcapital // OK
avgK_g1 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgK_g2 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgK_g3 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgK_g4 = (double *) calloc((ifulldimETsimul), sizeof(double));

// avgBANKRUPT // OK
avgBANKRUPT_g1 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgBANKRUPT_g2 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgBANKRUPT_g3 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgBANKRUPT_g4 = (double *) calloc((ifulldimETsimul), sizeof(double));

// avg MPLDL // OK
avgMPLDL_g1 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgMPLDL_g2 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgMPLDL_g3 = (double *) calloc((ifulldimETsimul), sizeof(double)); avgMPLDL_g4 = (double *) calloc((ifulldimETsimul), sizeof(double));

// Growth rate // OK
growthK_group1 = (double *) calloc((ifulldimETsimul), sizeof(double)); growthK_group2 = (double *) calloc((ifulldimETsimul), sizeof(double)); growthK_group3 = (double *) calloc((ifulldimETsimul), sizeof(double)); growthK_group4 = (double *) calloc((ifulldimETsimul), sizeof(double));

// Growth rate // OK
growthY_group1 = (double *) calloc((ifulldimETsimul), sizeof(double)); growthY_group2 = (double *) calloc((ifulldimETsimul), sizeof(double)); growthY_group3 = (double *) calloc((ifulldimETsimul), sizeof(double)); growthY_group4  = (double *) calloc((ifulldimETsimul), sizeof(double));
            
// SURVIVAL RATE // OK
surv_group1 = (double *) calloc((ifulldimETsimul), sizeof(double));
surv_group2 = (double *) calloc((ifulldimETsimul), sizeof(double));
surv_group3 = (double *) calloc((ifulldimETsimul), sizeof(double));
surv_group4 = (double *) calloc((ifulldimETsimul), sizeof(double));



// initialize integer //
int tt, ii, yy, ee, igrid, tgrid, zgrid, k, e, y, t, time;



#if SEA == 1
// VECTOR OF MUENDO FOR DRI //
double *MUENDOvecB, *MUENDOvecR, *MUENDOvecEhat;

MUENDOvecR      = (double *) calloc((ifulldimEE), sizeof(double));
MUENDOvecB      = (double *) calloc((ifulldimEE), sizeof(double));
MUENDOvecEhat   = (double *) calloc((ifulldimEE), sizeof(double));

for(igrid=0;igrid<maxgrid;igrid++){
    for(tgrid = 0.0; tgrid < maxtheta; tgrid ++){
        for(zgrid = 0.0; zgrid < maxfirmtype; zgrid++){
            for(e = 0; e < maxfirmtype; e++){
                
                MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)] = mupar*(bstarcostfun(profitRfun(igrid,e,tgrid,collateralJNEi[inxE(igrid,tgrid,zgrid)],rateJNEi[inxE(igrid,tgrid,zgrid)]),tgrid))/((1.0-taxrate)*rhostar*wstar*prod[tgrid]);
                
                MUENDOvecB[inxEE(igrid,tgrid,zgrid,e)] = mupar*(bstarcostfun(0.0,tgrid)/((1.0-taxrate)*rhostar*wstar*prod[tgrid]));
                
                MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)] = mupar*(bstarcostfun(profitEhatfun(igrid,e,tgrid,collateralJEi[inxE(igrid,tgrid,zgrid)]),tgrid))/((1.0-taxrate)*rhostar*wstar*prod[tgrid]);
                
                if(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)] < 0 || MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)] > mupar + 0.00000001){printf("mistake MUENDOvecEhat %f ", mupar*(bstarcostfun(profitEhatfun(igrid,zgrid,tgrid,collateralJEi[inxE(igrid,tgrid,zgrid)]),tgrid))/((1.0-taxrate)*rhostar*wstar*prod[tgrid]));getchar();}
                if(MUENDOvecB[inxEE(igrid,tgrid,zgrid,e)] < 0 || MUENDOvecB[inxEE(igrid,tgrid,zgrid,e)] > mupar + 0.00000001 ){printf("mistake MUENDOvecB %f", mupar);getchar();}
                if(MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)] < 0 || MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)] > mupar + 0.00000001){printf("mistake MUENDOvecR %f", mupar*(bstarcostfun(profitRfun(igrid,zgrid,tgrid,collateralJNEi[inxE(igrid,tgrid,zgrid)],rateJNEi[inxE(igrid,tgrid,zgrid)]),tgrid))/((1.0-taxrate)*rhostar*wstar*prod[tgrid]));getchar();}
            }
            
        }
    }
}

#endif


// INITIALIZE FIRST DISTRIBUTIONS //
for(ii = 0; ii < maxgrid; ii++)
{
    for(tt = 0; tt < maxtheta; tt++)
    {
        for(ee = 0; ee < maxfirmtype; ee++)
        {
            #if SEA == 1
            distEgroup1JNEitemp[inxETsimul(ii,tt,ee,0)] = newJNE1[inxE(ii,tt,ee)];
            distEgroup1JEitemp[inxETsimul(ii,tt,ee,0)] = newJE1[inxE(ii,tt,ee)];
            
            distEgroup2JNEitemp[inxETsimul(ii,tt,ee,0)] = newJNE2[inxE(ii,tt,ee)];
            distEgroup2JEitemp[inxETsimul(ii,tt,ee,0)] = newJE2[inxE(ii,tt,ee)];
            #endif
            
            #if NOPOLICY == 1
            distEgroup1JNEtemp[inxETsimul(ii,tt,ee,0)] = newJNE1[inxE(ii,tt,ee)];
            distEgroup1JEtemp[inxETsimul(ii,tt,ee,0)] = newJE1[inxE(ii,tt,ee)];
            
            distEgroup2JNEtemp[inxETsimul(ii,tt,ee,0)] = newJNE2[inxE(ii,tt,ee)];
            distEgroup2JEtemp[inxETsimul(ii,tt,ee,0)] = newJE2[inxE(ii,tt,ee)];
            #endif
    
            distEgroup3JNEtemp[inxETsimul(ii,tt,ee,0)] = newJNE3[inxE(ii,tt,ee)];
            distEgroup3JEtemp[inxETsimul(ii,tt,ee,0)] = newJE3[inxE(ii,tt,ee)];

//            distEgroup4JNEtemp[inxETsimul(ii,tt,ee,0)] = newJNE4[inxE(ii,tt,ee)];
//            distEgroup4JEtemp[inxETsimul(ii,tt,ee,0)] = newJE4[inxE(ii,tt,ee)];
        }
    }
}



double distverif1, distverif2;

distverif1 = 0.0;
distverif2 = 0.0;
double distverif3 = 0.0;
// COMPUTE STARTING DISTRIBUTIONS //
for(ii = 0; ii < maxgrid; ii++)
{
    for(tt = 0; tt < maxtheta; tt++)
    {
        for(ee = 0; ee < maxfirmtype; ee++)
        {
            #if SEA == 1
            dist_group1[inxETsimul(ii,tt,ee,0)] = distEgroup1JNEitemp[inxETsimul(ii, tt, ee, 0)]+distEgroup1JEitemp[inxETsimul(ii, tt, ee, 0)];
            dist_group2[inxETsimul(ii,tt,ee,0)] = (distEgroup2JNEitemp[inxETsimul(ii, tt, ee, 0)]+distEgroup2JEitemp[inxETsimul(ii, tt, ee, 0)]);
            #endif
            
            #if NOPOLICY == 1
            dist_group1[inxETsimul(ii,tt,ee,0)] = distEgroup1JNEtemp[inxETsimul(ii, tt, ee, 0)]+distEgroup1JEtemp[inxETsimul(ii, tt, ee, 0)];
            dist_group2[inxETsimul(ii,tt,ee,0)] = (distEgroup2JNEtemp[inxETsimul(ii, tt, ee, 0)]+distEgroup2JEtemp[inxETsimul(ii, tt, ee, 0)]);
            #endif
            
            dist_group3[inxETsimul(ii,tt,ee,0)] = (distEgroup3JEtemp[inxETsimul(ii, tt, ee, 0)]+distEgroup3JNEtemp[inxETsimul(ii,tt,ee, 0)]);
            dist_group4[inxETsimul(ii,tt,ee,0)] = (distEgroup4JEtemp[inxETsimul(ii, tt, ee, 0)]+distEgroup4JNEtemp[inxETsimul(ii,tt,ee, 0)]);


            distverif3 += dist_group3[inxETsimul(ii,tt,ee,0)];
            distverif2 += dist_group2[inxETsimul(ii,tt,ee,0)];
            distverif1 += dist_group1[inxETsimul(ii,tt,ee,0)];
        }
    }
}

printf("VERIF = %f %f %f %f %f %20.15f %20.15f %20.15f\n", distEgroup2JNEtemp[30], distEgroup2JNEtemp[30], distEgroup2JNEtemp[30], distEgroup2JEtemp[30], distEgroup2JEtemp[30], distverif1, distverif2,distverif3); // getchar();


for(t = 0; t < maxtimesimul; t++) {
    dist1[t] = 0.0;
    dist2[t] = 0.0;
    dist3[t] = 0.0;
    dist4[t] = 0.0;
    
    for(tt = 0; tt < maxtheta; tt++) {
            dist1T[tt*maxtimesimul + t] = 0.0;
            dist2T[tt*maxtimesimul + t] = 0.0;
            dist3T[tt*maxtimesimul + t] = 0.0;
            dist4T[tt*maxtimesimul + t] = 0.0;
    }
}

///////////////////////////////////////////////////////////////////////////
//////////////////////// SIMULATION FOR EACH POINT ////////////////////////
///////////////////////////////////////////////////////////////////////////

    
// UPPER LOOP, FOR EACH GRID POINT HAVE TO FIND THE AVERAGE OF EVERYTHING
for(ii = 0; ii < maxgrid; ii++)
{
    for(tt = 0; tt < maxtheta; tt++)
    {
        for(ee = 0; ee < maxfirmtype; ee++)
        {
            
            if((newJNE1[inxE(ii,tt,ee)] + newJE1[inxE(ii,tt,ee)] + newJNE2[inxE(ii,tt,ee)] + newJE2[inxE(ii,tt,ee)] + newJNE3[inxE(ii,tt,ee)] + newJE3[inxE(ii,tt,ee)]) > 0.0){
            
            
            // reinitialize vectors //
            bascule(zeroE, distEgroup1JNEi, ifulldimETsimul);
            bascule(zeroE, distEgroup1JNE, ifulldimETsimul);
            bascule(zeroE, distEgroup1JEi, ifulldimETsimul);
            bascule(zeroE, distEgroup1JE, ifulldimETsimul);

            bascule(zeroE, distEgroup2JNEi, ifulldimETsimul);
            bascule(zeroE, distEgroup2JNE, ifulldimETsimul);
            bascule(zeroE, distEgroup2JEi, ifulldimETsimul);
            bascule(zeroE, distEgroup2JE, ifulldimETsimul);
            
            bascule(zeroE, distEgroup3JNE, ifulldimETsimul);
            bascule(zeroE, distEgroup3JE, ifulldimETsimul);

            bascule(zeroE, distEgroup4JNE, ifulldimETsimul);
            bascule(zeroE, distEgroup4JE, ifulldimETsimul);
            
            
            // ALLOCATE VECTOR //
            // allocate vectors to compute the behavior at another grid point
            distEgroup1JNEi[inxETsimul(ii,tt,ee,0)] = distEgroup1JNEitemp[inxETsimul(ii,tt,ee,0)];
            distEgroup1JNE[inxETsimul(ii,tt,ee,0)] = distEgroup1JNEtemp[inxETsimul(ii,tt,ee,0)];
            distEgroup1JEi[inxETsimul(ii,tt,ee,0)] = distEgroup1JEitemp[inxETsimul(ii,tt,ee,0)];
            distEgroup1JE[inxETsimul(ii,tt,ee,0)] = distEgroup1JEtemp[inxETsimul(ii,tt,ee,0)];
            
            distEgroup2JNEi[inxETsimul(ii,tt,ee,0)] = distEgroup2JNEitemp[inxETsimul(ii,tt,ee,0)];
            distEgroup2JNE[inxETsimul(ii,tt,ee,0)] = distEgroup2JNEtemp[inxETsimul(ii,tt,ee,0)];
            distEgroup2JEi[inxETsimul(ii,tt,ee,0)] = distEgroup2JEitemp[inxETsimul(ii,tt,ee,0)];
            distEgroup2JE[inxETsimul(ii,tt,ee,0)] = distEgroup2JEtemp[inxETsimul(ii,tt,ee,0)];
            
            distEgroup3JNE[inxETsimul(ii,tt,ee,0)] = distEgroup3JNEtemp[inxETsimul(ii,tt,ee,0)];
            distEgroup3JE[inxETsimul(ii,tt,ee,0)] = distEgroup3JEtemp[inxETsimul(ii,tt,ee,0)];
            
            distEgroup4JNE[inxETsimul(ii,tt,ee,0)] = distEgroup4JNEtemp[inxETsimul(ii,tt,ee,0)];
            distEgroup4JE[inxETsimul(ii,tt,ee,0)] = distEgroup4JEtemp[inxETsimul(ii,tt,ee,0)];
            
            
            // STARTING HERE THE BIG GUYS //
            // start by a upper loop on time //
            for(t = 0; t < maxtimesimul; t++)
            {
            
   
//            if(ii == 20 && tt == 1 && ee == 1) {
//             printf("%f  ", 10000*(distEgroup3JNE[inxETsimul(ii,tt,ee,t-1)] + distEgroup3JE[inxETsimul(ii,tt,ee,t-1)])); getchar();
//            }

            for(igrid = 0; igrid < maxgrid; igrid++)
            {
            for(tgrid = 0; tgrid < maxtheta; tgrid++)
            {
            for(k = 0; k < maxtheta; k++)
            {
            for(zgrid = 0; zgrid < maxfirmtype; zgrid++)
            {
            for(e = 0; e < maxfirmtype; e++)
            {
            
            if(t>0){
            
                #if SEA == 1
                
                /*****************************************/
                /*** ENTREPRENEUR INSURED NOT-EXCLUDED ***/
                /*****************************************/


                // if the entrepreneur repay
                if(defaultJNEi[inxEE(igrid, tgrid, zgrid, e)] == 0) {
                    
                    // TRANSITION ENT TO WORKER NOT EXCLUDED (receive job offer)
                    for(y=0; y < maxyprod; y++) {

                        // does not lost his rights // residual who stay entrepreneur
                        expprob = invmyprod[y]*mprod[tgrid][k]*(1.0 - MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)])*mzprob[zgrid][e]*(piW( searchJNEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                        distEgroup1JNEi[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup1JNEi[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                        
                        distEgroup2JNEi[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup2JNEi[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                        

                        // lost his rights // residual who stay entrepreneur
                        expprob = invmyprod[y]*mprod[tgrid][k]*(MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)])*mzprob[zgrid][e]*(piW( searchJNEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                        distEgroup1JNE[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                        
                        distEgroup2JNE[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
            
                    }
                

                    // residual who stay entrepreneur
                    expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - piW( searchJNEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JNEi[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUiNE[inx(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup1JNEi[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUiNE[inx(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - piW( searchJNEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup2JNEi[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUiNE[inx(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup2JNEi[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUiNE[inx(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    
                    // residual who stay entrepreneur
                    expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecR[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - piW( searchJNEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JNE[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    distEgroup2JNE[inxETsimul(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJNEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    
                } // end condition




                /*************************************/
                /*** ENTREPRENEUR INSURED EXCLUDED ***/
                /*************************************/

                for(y=0; y < maxyprod; y++) {

                    
                    // and residual who stay entrepreneur
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - phipar)*(piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup1JEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - phipar)*(piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup2JEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup2JEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    
                    // and residual who stay entrepreneur
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(phipar)*(piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JNEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup1JNEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup2JNEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup2JNEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    


                    // and residual who stay entrepreneur
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - phipar)*(piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup1JE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup2JE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup2JE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    

                    
                    // and residual who stay entrepreneur
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(phipar)*(piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JNE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup2JNE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    
                }


                    // and residual who star entrepreneur
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - phipar)*(1.0 - piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup1JEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUiE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup1JEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUiE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup2JEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUiE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup2JEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUiE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                
                
                // and residual who star entrepreneur
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(phipar)*(1.0 - piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup1JNEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUiNE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup1JNEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUiNE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup2JNEi[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUiNE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup2JNEi[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNEi[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUiNE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));




                // TRANSITION ENT TO UNEMPLOYED LT EXCLUDED  // lose his rights
                // TRANSITION ENT TO UNINSURED UNEMPLOYED EXCLUDED (receive no job offer, not receive phipar)

                // and residual who star entrepreneur
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - phipar)*(1.0 - piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup1JE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup1JE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(1.0 - phipar)*(1.0 - piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup2JE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup2JE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                
                
                // and residual who star entrepreneur
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(phipar)*(1.0 - piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup1JNE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(MUENDOvecEhat[inxEE(igrid,tgrid,zgrid,e)])*(phipar)*(1.0 - piW( searchJEi_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup2JNE[inxETsimul(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJEi[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJEi[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJEi[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                
                
                #endif
                
                
                

                /*********************************/
                /*** ENTREPRENEUR NOT-EXCLUDED ***/
                /*********************************/
                
                // if the entrepreneur repay
                if(defaultJNE[inxEE(igrid, tgrid, zgrid, e)] == 0) {
                    
                    // TRANSITION ENT TO WORKER NOT EXCLUDED (receive job offer)
                    for(y=0; y < maxyprod; y++) {
                    
                        // residual who stay entrepreneur
                        expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(piW( searchJNE_W[inxEE(igrid,tgrid,zgrid, e)]));
                        distEgroup1JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                        
                        distEgroup2JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                        
                        distEgroup3JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup3JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                        
                        distEgroup4JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                        distEgroup4JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));

                    }
                    
                    // residual who stay entrepreneur
                    expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - piW( searchJNE_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    distEgroup2JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    distEgroup3JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup3JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                    
                    distEgroup4JNE[inxETsimul(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJNE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                    distEgroup4JNE[inxETsimul(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJNE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJNE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
    
                    
                } // end condition
                

                 

                
                /*****************************/
                /*** ENTREPRENEUR EXCLUDED ***/
                /*****************************/
                
                for(y=0; y < maxyprod; y++) {

                    // and residual who stay entrepreneur
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - phipar)*(piW( searchJE_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup1JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup2JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup2JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup3JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup3JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup4JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup4JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    // and residual who stay entrepreneur
                    expprob = invmyprod[y]*mprod[tgrid][k]*mzprob[zgrid][e]*(phipar)*(piW( searchJE_W[inxEE(igrid,tgrid,zgrid, e)]));
                    distEgroup1JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup2JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup3JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup3JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                    
                    distEgroup4JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e, t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] > (valueWWNE[inxW(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, y)])))?(1.0):(0.0));
                    distEgroup4JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e, t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] > (valueWWNE[inxW(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, y)])))?(1.0):(0.0));
                }
                

                // and residual who star entrepreneur
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(1.0 - phipar)*(1.0 - piW( searchJE_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup1JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup1JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup2JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup2JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup3JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup3JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup4JE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup4JE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));


                // and residual who star entrepreneur
                expprob = mprod[tgrid][k]*mzprob[zgrid][e]*(phipar)*(1.0 - piW( searchJE_W[inxEE(igrid,tgrid,zgrid, e)]));
                distEgroup1JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup1JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup2JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup2JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup3JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup3JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
                distEgroup4JNE[inxETsimul(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e,t)] += expprob*(1.0 - saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k, e)] >= (valueUnNE[inx(isaveJE[inxEE(igrid, tgrid, zgrid, e)], k)])))?(1.0):(0.0));
                distEgroup4JNE[inxETsimul(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e,t)] += expprob*(saveresJE[inxEE(igrid, tgrid, zgrid, e)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]*(((valueJNE[inxE(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k, e)] >= (valueUnNE[inx(min((maxgrid-1),(isaveJE[inxEE(igrid, tgrid, zgrid, e)]+1)), k)])))?(1.0):(0.0));
                
            } // end e
            } // end zgrid
            } // end k
            } // end tgrid
            } // end igrid
            } // end if t > 0

            // avg wealth level //
            avgW_group1 = 0.0;
            avgW_group2 = 0.0;
            avgW_group3 = 0.0;
            avgW_group4 = 0.0;
            
            // Total dist //
            distgroup1 = 0.0;
            distgroup2 = 0.0;
            distgroup3 = 0.0;
            distgroup4 = 0.0;
            
            // Average skill //
            avgT_group1 = 0.0;
            avgT_group2 = 0.0;
            avgT_group3 = 0.0;
            avgT_group4 = 0.0;
                
            // Average prod // OK
            avgY_group1 = 0.0;
            avgY_group2 = 0.0;
            avgY_group3 = 0.0;
            avgY_group4 = 0.0;
                
            // avgcapital // OK
            avgK_group1 = 0.0;
            avgK_group2 = 0.0;
            avgK_group3 = 0.0;
            avgK_group4 = 0.0;
                
            // avgBANKRUPT // OK
            avgBANKRUPT_group1 = 0.0;
            avgBANKRUPT_group2 = 0.0;
            avgBANKRUPT_group3 = 0.0;
            avgBANKRUPT_group4 = 0.0;
            
            // avg MPLDL // OK
            avgMPLDL_group1 = 0.0;
            avgMPLDL_group2 = 0.0;
            avgMPLDL_group3 = 0.0;
            avgMPLDL_group4 = 0.0;
                
            // Average prod // OK
            avgY_group1_bef = 0.0;
            avgY_group2_bef = 0.0;
            avgY_group3_bef = 0.0;
            avgY_group4_bef = 0.0;
                
            // avgcapital // OK
            avgK_group1_bef = 0.0;
            avgK_group2_bef = 0.0;
            avgK_group3_bef = 0.0;
            avgK_group4_bef = 0.0;
            
            
            
            for(igrid = 0; igrid < maxgrid; igrid++){
            for(tgrid = 0; tgrid < maxtheta; tgrid++){
            for(zgrid = 0; zgrid < maxfirmtype; zgrid++){
            
            
            // HERE COMPUTE INDICATOR FOR EACH GROUP FOR PERIOD T, WITH INITIAL CONDITION IGRID, TGRID, ZGRID //
            
            // avg wealth level //
            avgW_group1 += grid[igrid]*(distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgW_group2 += grid[igrid]*(distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgW_group3 += grid[igrid]*(distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgW_group4 += grid[igrid]*(distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
            
                
           // avg bankruptcy rate //
           for(e = 0; e < maxfirmtype; e++)
           {
           #if SEA == 1
           avgBANKRUPT_group1 += mzprob[zgrid][e]*(defaultJNEi[inxEE(igrid, tgrid, zgrid, e)]*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+defaultJNE[inxEE(igrid, tgrid, zgrid, e)]*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
           avgBANKRUPT_group2 += mzprob[zgrid][e]*(defaultJNEi[inxEE(igrid, tgrid, zgrid, e)]*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+defaultJNE[inxEE(igrid, tgrid, zgrid, e)]*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
           #endif
           #if NOPOLICY == 1 
           avgBANKRUPT_group1 += mzprob[zgrid][e]*(defaultJNE[inxEE(igrid, tgrid, zgrid, e)]*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
           avgBANKRUPT_group2 += mzprob[zgrid][e]*(defaultJNE[inxEE(igrid, tgrid, zgrid, e)]*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
           #endif
           avgBANKRUPT_group3 += mzprob[zgrid][e]*(defaultJNE[inxEE(igrid, tgrid, zgrid, e)]*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
           avgBANKRUPT_group4 += mzprob[zgrid][e]*(defaultJNE[inxEE(igrid, tgrid, zgrid, e)]*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
           }
                
                
            // Total dist //
            distgroup1 += (distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            distgroup2 += (distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            distgroup3 += (distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
            distgroup4 += (distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
            
            
            // Average skill //
            avgT_group1 += gmapping[tgrid]*(distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgT_group2 += gmapping[tgrid]*(distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgT_group3 += gmapping[tgrid]*(distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgT_group4 += gmapping[tgrid]*(distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
                
                
            // Average prod // OK
            #if SEA == 1
            avgY_group1 += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNEi[inxE(igrid,tgrid,ee)]))*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJEi[inxE(igrid,tgrid,ee)]))*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t)];
            
            avgY_group2 += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNEi[inxE(igrid,tgrid,ee)]))*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJEi[inxE(igrid,tgrid,ee)]))*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t)];
            #endif
            
            #if NOPOLICY == 1
            avgY_group1 += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            
            avgY_group2 += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            #endif
            
            avgY_group3 += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            
            avgY_group4 += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)];
                
                
            // avgcapital // OK
            #if SEA == 1
            avgK_group1 += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNEi[inxE(igrid, tgrid, zgrid)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJEi[inxE(igrid, tgrid, zgrid)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t)];
            
            avgK_group2 += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNEi[inxE(igrid, tgrid, zgrid)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJEi[inxE(igrid, tgrid, zgrid)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t)];
            #endif
            
            #if NOPOLICY == 1
            avgK_group1 += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            
            avgK_group2 += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            #endif
            
            avgK_group3 += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            
            avgK_group4 += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)];
            
            
            // avg MPLDL // OK
            avgMPLDL_group1 += ((1-alphapar)*TFPpar*pow((aggcapital)/(agglabor),alphapar)*prod[tgrid])*(distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgMPLDL_group2 += ((1-alphapar)*TFPpar*pow((aggcapital)/(agglabor),alphapar)*prod[tgrid])*(distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgMPLDL_group3 += ((1-alphapar)*TFPpar*pow((aggcapital)/(agglabor),alphapar)*prod[tgrid])*(distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
            avgMPLDL_group4 += ((1-alphapar)*TFPpar*pow((aggcapital)/(agglabor),alphapar)*prod[tgrid])*(distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t)]+distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t)]);
                
            if(t>0){
            // Average prod bef// OK
            #if NOPOLICY == 1
            avgY_group1_bef += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            
            avgY_group2_bef += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            #endif
            
            #if SEA == 1
            avgY_group1_bef += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNEi[inxE(igrid,tgrid,ee)]))*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJEi[inxE(igrid,tgrid,ee)]))*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)];
            
            avgY_group2_bef += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNEi[inxE(igrid,tgrid,ee)]))*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJEi[inxE(igrid,tgrid,ee)]))*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)];
            #endif
            
            avgY_group3_bef += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            
            avgY_group4_bef += (prodE(statez[zgrid],tgrid,collateralJE[inxE(igrid,tgrid,ee)]))*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(prodE(statez[zgrid],tgrid,collateralJNE[inxE(igrid,tgrid,ee)]))*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
                
                
            // avgcapital bef // OK
            #if SEA == 1
            avgK_group1_bef += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNEi[inxE(igrid, tgrid, zgrid)])*distEgroup1JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJEi[inxE(igrid, tgrid, zgrid)])*distEgroup1JEi[inxETsimul(igrid, tgrid, zgrid, t-1)];
            
            avgK_group2_bef += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNEi[inxE(igrid, tgrid, zgrid)])*distEgroup2JNEi[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJEi[inxE(igrid, tgrid, zgrid)])*distEgroup2JEi[inxETsimul(igrid, tgrid, zgrid, t-1)];
            #endif
            
            #if NOPOLICY == 1 
            avgK_group1_bef += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup1JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup1JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            
            avgK_group2_bef += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup2JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup2JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            #endif
            
            avgK_group3_bef += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup3JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup3JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            
            avgK_group4_bef += (collateralJE[inxE(igrid, tgrid, zgrid)])*distEgroup4JE[inxETsimul(igrid, tgrid, zgrid, t-1)]+(collateralJNE[inxE(igrid, tgrid, zgrid)])*distEgroup4JNE[inxETsimul(igrid, tgrid, zgrid, t-1)];
            }

            
            } // end zgrid
            } // end tgrid
            } // end igrid
            
            
            // DIST //
            dist_group1[inxETsimul(ii,tt,ee,t)] = distgroup1;
            dist_group2[inxETsimul(ii,tt,ee,t)] = distgroup2;
            dist_group3[inxETsimul(ii,tt,ee,t)] = distgroup3;
            dist_group4[inxETsimul(ii,tt,ee,t)] = distgroup4;
            
                
            // AVG W //
            avgW_g1[inxETsimul(ii,tt,ee,t)] = avgW_group1;
            avgW_g2[inxETsimul(ii,tt,ee,t)] = avgW_group2;
            avgW_g3[inxETsimul(ii,tt,ee,t)] = avgW_group3;
            avgW_g4[inxETsimul(ii,tt,ee,t)] = avgW_group4;

            // AVG T //
            avgT_g1[inxETsimul(ii,tt,ee,t)] = avgT_group1;
            avgT_g2[inxETsimul(ii,tt,ee,t)] = avgT_group2;
            avgT_g3[inxETsimul(ii,tt,ee,t)] = avgT_group3;
            avgT_g4[inxETsimul(ii,tt,ee,t)] = avgT_group4;
                
            // AVG Y //
            avgY_g1[inxETsimul(ii,tt,ee,t)] = avgY_group1;
            avgY_g2[inxETsimul(ii,tt,ee,t)] = avgY_group2;
            avgY_g3[inxETsimul(ii,tt,ee,t)] = avgY_group3;
            avgY_g4[inxETsimul(ii,tt,ee,t)] = avgY_group4;
               
            // AVG K //
            avgK_g1[inxETsimul(ii,tt,ee,t)] = avgK_group1;
            avgK_g2[inxETsimul(ii,tt,ee,t)] = avgK_group2;
            avgK_g3[inxETsimul(ii,tt,ee,t)] = avgK_group3;
            avgK_g4[inxETsimul(ii,tt,ee,t)] = avgK_group4;
                
            // AVG BANKRUPT //
            avgBANKRUPT_g1[inxETsimul(ii,tt,ee,t)] = avgBANKRUPT_group1;
            avgBANKRUPT_g2[inxETsimul(ii,tt,ee,t)] = avgBANKRUPT_group2;
            avgBANKRUPT_g3[inxETsimul(ii,tt,ee,t)] = avgBANKRUPT_group3;
            avgBANKRUPT_g4[inxETsimul(ii,tt,ee,t)] = avgBANKRUPT_group4;

            // AVG MPLDL //
            avgMPLDL_g1[inxETsimul(ii,tt,ee,t)] = avgMPLDL_group1;
            avgMPLDL_g2[inxETsimul(ii,tt,ee,t)] = avgMPLDL_group2;
            avgMPLDL_g3[inxETsimul(ii,tt,ee,t)] = avgMPLDL_group3;
            avgMPLDL_g4[inxETsimul(ii,tt,ee,t)] = avgMPLDL_group4;
              
            if(t>0){
            // Growth rate Y // OK
            if(avgY_group1_bef > 0.0 && dist_group1[inxETsimul(ii,tt,ee,t-1)] > 0.0 && dist_group1[inxETsimul(ii,tt,ee,t)] > 0.0){growthY_group1[inxETsimul(ii,tt,ee,t)] = (avgY_group1/dist_group1[inxETsimul(ii,tt,ee,t)] - avgY_group1_bef/dist_group1[inxETsimul(ii,tt,ee,t-1)])/(avgY_group1_bef/dist_group1[inxETsimul(ii,tt,ee,t-1)]);}
            if(avgY_group2_bef > 0.0 && dist_group2[inxETsimul(ii,tt,ee,t-1)] > 0.0 && dist_group2[inxETsimul(ii,tt,ee,t)] > 0.0){growthY_group2[inxETsimul(ii,tt,ee,t)] = (avgY_group2/dist_group2[inxETsimul(ii,tt,ee,t)] - avgY_group2_bef/dist_group2[inxETsimul(ii,tt,ee,t-1)])/(avgY_group2_bef/dist_group2[inxETsimul(ii,tt,ee,t-1)]);}
            if(avgY_group3_bef > 0.0 && dist_group3[inxETsimul(ii,tt,ee,t-1)] > 0.0  && dist_group3[inxETsimul(ii,tt,ee,t)] > 0.0){growthY_group3[inxETsimul(ii,tt,ee,t)] = (avgY_group3/dist_group3[inxETsimul(ii,tt,ee,t)] - avgY_group3_bef/dist_group3[inxETsimul(ii,tt,ee,t-1)])/(avgY_group3_bef/dist_group3[inxETsimul(ii,tt,ee,t-1)]);}
            if(avgY_group4_bef > 0.0 && dist_group4[inxETsimul(ii,tt,ee,t-1)] > 0.0  && dist_group4[inxETsimul(ii,tt,ee,t)] > 0.0){growthY_group4[inxETsimul(ii,tt,ee,t)] = (avgY_group4/dist_group4[inxETsimul(ii,tt,ee,t)] - avgY_group4_bef/dist_group4[inxETsimul(ii,tt,ee,t-1)])/(avgY_group4_bef/dist_group4[inxETsimul(ii,tt,ee,t-1)]);}


            // Growth rate K // OK
            if(avgK_group1_bef > 0.0 && dist_group1[inxETsimul(ii,tt,ee,t-1)] > 0.0 && dist_group1[inxETsimul(ii,tt,ee,t)] > 0.0){growthK_group1[inxETsimul(ii,tt,ee,t)] = (avgK_group1/dist_group1[inxETsimul(ii,tt,ee,t)] - avgK_group1_bef/dist_group1[inxETsimul(ii,tt,ee,t-1)])/(avgK_group1_bef/dist_group1[inxETsimul(ii,tt,ee,t-1)]);}
            if(avgK_group2_bef > 0.0 && dist_group2[inxETsimul(ii,tt,ee,t-1)] > 0.0 && dist_group2[inxETsimul(ii,tt,ee,t)] >0.0){growthK_group2[inxETsimul(ii,tt,ee,t)] = (avgK_group2/dist_group2[inxETsimul(ii,tt,ee,t)] - avgK_group2_bef/dist_group2[inxETsimul(ii,tt,ee,t-1)])/(avgK_group2_bef/dist_group2[inxETsimul(ii,tt,ee,t-1)]);}
            if(avgK_group3_bef > 0.0 && dist_group3[inxETsimul(ii,tt,ee,t-1)] > 0.0 && dist_group3[inxETsimul(ii,tt,ee,t)] > 0.0){growthK_group3[inxETsimul(ii,tt,ee,t)] = (avgK_group3/dist_group3[inxETsimul(ii,tt,ee,t)] - avgK_group3_bef/dist_group3[inxETsimul(ii,tt,ee,t-1)])/(avgK_group3_bef/dist_group3[inxETsimul(ii,tt,ee,t-1)]);}
            if(avgK_group4_bef > 0.0 && dist_group4[inxETsimul(ii,tt,ee,t-1)] > 0.0 && dist_group4[inxETsimul(ii,tt,ee,t)] > 0.0){growthK_group4[inxETsimul(ii,tt,ee,t)] = (avgK_group4/dist_group4[inxETsimul(ii,tt,ee,t)] - avgK_group4_bef/dist_group4[inxETsimul(ii,tt,ee,t-1)])/(avgK_group4_bef/dist_group4[inxETsimul(ii,tt,ee,t-1)]);}
            

            // survival rate //
            if(dist_group1[inxETsimul(ii,tt,ee,0)] > 0.0){surv_group1[inxETsimul(ii,tt,ee,t)] = distgroup1/dist_group1[inxETsimul(ii,tt,ee,0)];}
            if(dist_group2[inxETsimul(ii,tt,ee,0)] > 0.0){surv_group2[inxETsimul(ii,tt,ee,t)] = distgroup2/dist_group2[inxETsimul(ii,tt,ee,0)];}
            if(dist_group3[inxETsimul(ii,tt,ee,0)] > 0.0){surv_group3[inxETsimul(ii,tt,ee,t)] = distgroup3/dist_group3[inxETsimul(ii,tt,ee,0)];}
            if(dist_group4[inxETsimul(ii,tt,ee,0)] > 0.0){surv_group4[inxETsimul(ii,tt,ee,t)] = distgroup4/dist_group4[inxETsimul(ii,tt,ee,0)];}
            }
            
            // DISTRIBUTION //
            dist1[t] += dist_group1[inxETsimul(ii,tt,ee,t)];
            dist2[t] += dist_group2[inxETsimul(ii,tt,ee,t)];
            dist3[t] += dist_group3[inxETsimul(ii,tt,ee,t)];
            dist4[t] += dist_group4[inxETsimul(ii,tt,ee,t)];
            
            // DISTRIBUTION THETA//
            dist1T[tt*maxtimesimul + t] += distgroup1;
            dist2T[tt*maxtimesimul + t] += distgroup2;
            dist3T[tt*maxtimesimul + t] += distgroup3;
            dist4T[tt*maxtimesimul + t] += distgroup4;
                
            } // end time
            
            }
            
            printf("\r                                              ");
            printf("\rComputation: %d over %d\n", inxE(ii,tt,ee), ifulldimE);
        } // end ee
    } // end tt
} // end ii








//////////////////////////////////////////////////////////////////////////////
//////////////////////// COMPUTE AVERAGE CONTRIBUTION ////////////////////////
//////////////////////////////////////////////////////////////////////////////


double avgW1[maxtimesimul], avgW2[maxtimesimul], avgW3[maxtimesimul], avgW4[maxtimesimul],
avgT1[maxtimesimul], avgT2[maxtimesimul], avgT3[maxtimesimul], avgT4[maxtimesimul],
avgY1[maxtimesimul], avgY2[maxtimesimul], avgY3[maxtimesimul], avgY4[maxtimesimul], 
avgK1[maxtimesimul], avgK2[maxtimesimul], avgK3[maxtimesimul], avgK4[maxtimesimul], 
avgBANK1[maxtimesimul], avgBANK2[maxtimesimul], avgBANK3[maxtimesimul], avgBANK4[maxtimesimul],
avgMPLDL1[maxtimesimul], avgMPLDL2[maxtimesimul], avgMPLDL3[maxtimesimul], avgMPLDL4[maxtimesimul],
growthK1[maxtimesimul], growthK2[maxtimesimul], growthK3[maxtimesimul], growthK4[maxtimesimul],
growthY1[maxtimesimul], growthY2[maxtimesimul], growthY3[maxtimesimul], growthY4[maxtimesimul],
surv1_byT[maxtimesimul*maxtheta], surv2_byT[maxtimesimul*maxtheta], surv3_byT[maxtimesimul*maxtheta], surv4_byT[maxtimesimul*maxtheta],
surv1[maxtimesimul], surv2[maxtimesimul], surv3[maxtimesimul], surv4[maxtimesimul];


double avgtW1,  avgtW2,  avgtW3,  avgtW4,
avgtT1,  avgtT2,  avgtT3,  avgtT4, 
avgtY1,  avgtY2,  avgtY3,  avgtY4,  
avgtK1,  avgtK2,  avgtK3,  avgtK4,  
avgtBK1,  avgtBK2,  avgtBK3,  avgtBK4,
avgtMPLDL1,  avgtMPLDL2,  avgtMPLDL3,  avgtMPLDL4,
growthtK1,  growthtK2,  growthtK3,  growthtK4,
growthtY1,  growthtY2,  growthtY3,  growthtY4,
surv1t, surv2t, surv3t, surv4t;

for(t = 0; t < maxtimesimul; t++){
    avgW1[t] = 0.0; avgW2[t] = 0.0; avgW3[t] = 0.0; avgW4[t] = 0.0;
    avgT1[t] = 0.0; avgT2[t] = 0.0; avgT3[t] = 0.0; avgT4[t] = 0.0;
    avgY1[t] = 0.0; avgY2[t] = 0.0; avgY3[t] = 0.0; avgY4[t] = 0.0; 
    avgK1[t] = 0.0; avgK2[t] = 0.0; avgK3[t] = 0.0; avgK4[t] = 0.0;
    avgBANK1[t] = 0.0; avgBANK2[t] = 0.0; avgBANK3[t] = 0.0; avgBANK4[t] = 0.0;
    avgMPLDL1[t] = 0.0; avgMPLDL2[t] = 0.0; avgMPLDL3[t] = 0.0; avgMPLDL4[t] = 0.0;
    growthK1[t] = 0.0; growthK2[t] = 0.0; growthK3[t] = 0.0; growthK4[t] = 0.0;
    growthY1[t] = 0.0; growthY2[t] = 0.0; growthY3[t] = 0.0; growthY4[t] = 0.0;
    surv1[t] = 0.0; surv2[t] = 0.0; surv3[t] = 0.0; surv4[t] = 0.0;
    
    for(tt = 0; tt < maxtheta; tt++){
        surv1_byT[tt*maxtimesimul + t] = 0.0;
        surv2_byT[tt*maxtimesimul + t] = 0.0;
        surv3_byT[tt*maxtimesimul + t] = 0.0;
        surv4_byT[tt*maxtimesimul + t] = 0.0;
    }
}

avgtW1 = 0.0; avgtW2 = 0.0; avgtW3 = 0.0; avgtW4 = 0.0;
avgtT1 = 0.0; avgtT2 = 0.0; avgtT3 = 0.0; avgtT4 = 0.0;
avgtY1 = 0.0; avgtY2 = 0.0; avgtY3 = 0.0; avgtY4 = 0.0; 
avgtK1 = 0.0; avgtK2 = 0.0; avgtK3 = 0.0; avgtK4 = 0.0;
avgtBK1 = 0.0; avgtBK2 = 0.0; avgtBK3 = 0.0; avgtBK4 = 0.0;
surv1t = 0.0; surv2t = 0.0; surv3t = 0.0; surv4t = 0.0;
avgtMPLDL1 = 0.0; avgtMPLDL2 = 0.0; avgtMPLDL3 = 0.0; avgtMPLDL4 = 0.0; 
growthtK1 = 0.0; growthtK2 = 0.0; growthtK3 = 0.0; growthtK4 = 0.0;
growthtY1 = 0.0; growthtY2 = 0.0; growthtY3 = 0.0; growthtY4 = 0.0;
            
    
    


// COMPUTE AVERAGE PER YEAR AND OVER 10 YEARS //
for(ii = 0; ii < maxgrid; ii++){
    for(tt = 0; tt < maxtheta; tt++){
        for(ee = 0; ee < maxfirmtype; ee++){
            for(t = 0; t < maxtimesimul; t++){
    
            
            // AVG WEALTH //
            avgW1[t] += avgW_g1[inxETsimul(ii,tt,ee,t)]/dist1[t];
            avgW2[t] += avgW_g2[inxETsimul(ii,tt,ee,t)]/dist2[t];
            avgW3[t] += avgW_g3[inxETsimul(ii,tt,ee,t)]/dist3[t];
            avgW4[t] += avgW_g4[inxETsimul(ii,tt,ee,t)]/dist4[t];
            
            
            // AVG SKILL //
            avgT1[t] += avgT_g1[inxETsimul(ii,tt,ee,t)]/dist1[t];
            avgT2[t] += avgT_g2[inxETsimul(ii,tt,ee,t)]/dist2[t];
            avgT3[t] += avgT_g3[inxETsimul(ii,tt,ee,t)]/dist3[t];
            avgT4[t] += avgT_g4[inxETsimul(ii,tt,ee,t)]/dist4[t];
            
            
            //AVG Y //
            avgY1[t] += avgY_g1[inxETsimul(ii,tt,ee,t)]/dist1[t];
            avgY2[t] += avgY_g2[inxETsimul(ii,tt,ee,t)]/dist2[t];
            avgY3[t] += avgY_g3[inxETsimul(ii,tt,ee,t)]/dist3[t];
            avgY4[t] += avgY_g4[inxETsimul(ii,tt,ee,t)]/dist4[t];
            
            
            // AVG K //
            avgK1[t] += avgK_g1[inxETsimul(ii,tt,ee,t)]/dist1[t];
            avgK2[t] += avgK_g2[inxETsimul(ii,tt,ee,t)]/dist2[t];
            avgK3[t] += avgK_g3[inxETsimul(ii,tt,ee,t)]/dist3[t];
            avgK4[t] += avgK_g4[inxETsimul(ii,tt,ee,t)]/dist4[t];
            
            // AVG K //
            avgBANK1[t] += avgBANKRUPT_g1[inxETsimul(ii,tt,ee,t)]/dist1[t];
            avgBANK2[t] += avgBANKRUPT_g2[inxETsimul(ii,tt,ee,t)]/dist2[t];
            avgBANK3[t] += avgBANKRUPT_g3[inxETsimul(ii,tt,ee,t)]/dist3[t];
            avgBANK4[t] += avgBANKRUPT_g4[inxETsimul(ii,tt,ee,t)]/dist4[t];
                
            // AVG MPLDL //
            avgMPLDL1[t] += avgMPLDL_g1[inxETsimul(ii,tt,ee,t)]/dist1[t];
            avgMPLDL2[t] += avgMPLDL_g2[inxETsimul(ii,tt,ee,t)]/dist2[t];
            avgMPLDL3[t] += avgMPLDL_g3[inxETsimul(ii,tt,ee,t)]/dist3[t];
            avgMPLDL4[t] += avgMPLDL_g4[inxETsimul(ii,tt,ee,t)]/dist4[t];
            
            
            // AVG GROWTH K //
            growthK1[t] += (growthK_group1[inxETsimul(ii,tt,ee,t)]*dist_group1[inxETsimul(ii,tt,ee,t)])/dist1[t];
            growthK2[t] += (growthK_group2[inxETsimul(ii,tt,ee,t)]*dist_group2[inxETsimul(ii,tt,ee,t)])/dist2[t];
            growthK3[t] += (growthK_group3[inxETsimul(ii,tt,ee,t)]*dist_group3[inxETsimul(ii,tt,ee,t)])/dist3[t];
            growthK4[t] += (growthK_group4[inxETsimul(ii,tt,ee,t)]*dist_group4[inxETsimul(ii,tt,ee,t)])/dist4[t];
            
            if(t>0){
            // AVG GROWTH Y //
            growthY1[t] += (growthY_group1[inxETsimul(ii,tt,ee,t)]*dist_group1[inxETsimul(ii,tt,ee,t)])/dist1[t];
            growthY2[t] += (growthY_group2[inxETsimul(ii,tt,ee,t)]*dist_group2[inxETsimul(ii,tt,ee,t)])/dist2[t];
            growthY3[t] += (growthY_group3[inxETsimul(ii,tt,ee,t)]*dist_group3[inxETsimul(ii,tt,ee,t)])/dist3[t];
            growthY4[t] += (growthY_group4[inxETsimul(ii,tt,ee,t)]*dist_group4[inxETsimul(ii,tt,ee,t)])/dist4[t];
        
            
            // survival rate //
            surv1_byT[tt*maxtimesimul + t] += (surv_group1[inxETsimul(ii,tt,ee,t)]*dist_group1[inxETsimul(ii,tt,ee,t)])/dist1T[tt*maxtimesimul + t];
            surv2_byT[tt*maxtimesimul + t] += (surv_group2[inxETsimul(ii,tt,ee,t)]*dist_group2[inxETsimul(ii,tt,ee,t)])/dist2T[tt*maxtimesimul + t];
            surv3_byT[tt*maxtimesimul + t] += (surv_group3[inxETsimul(ii,tt,ee,t)]*dist_group3[inxETsimul(ii,tt,ee,t)])/dist3T[tt*maxtimesimul + t];
            surv4_byT[tt*maxtimesimul + t] += (surv_group4[inxETsimul(ii,tt,ee,t)]*dist_group4[inxETsimul(ii,tt,ee,t)])/dist4T[tt*maxtimesimul + t];
            
            surv1[t] += (surv_group1[inxETsimul(ii,tt,ee,t)]*dist_group1[inxETsimul(ii,tt,ee,t)])/dist1[t];
            surv2[t] += (surv_group2[inxETsimul(ii,tt,ee,t)]*dist_group2[inxETsimul(ii,tt,ee,t)])/dist2[t];
            surv3[t] += (surv_group3[inxETsimul(ii,tt,ee,t)]*dist_group3[inxETsimul(ii,tt,ee,t)])/dist3[t];
            surv4[t] += (surv_group4[inxETsimul(ii,tt,ee,t)]*dist_group4[inxETsimul(ii,tt,ee,t)])/dist4[t];
            }
        

            } // end t
        } // end ee
    } // end tt
} // end ii



distverif1 = 0.0;
distverif2 = 0.0;
// COMPUTE STARTING DISTRIBUTIONS //
for(ii = 0; ii < maxgrid; ii++)
{
    for(tt = 0; tt < maxtheta; tt++)
    {
        for(ee = 0; ee < maxfirmtype; ee++)
        {
            distverif1 += dist_group1[inxETsimul(ii,tt,ee,(maxtime-1))];
            distverif2 += dist_group2[inxETsimul(ii,tt,ee,(maxtime-1))];
        }
    }
}

printf("VERIF = %20.15f %20.15f\n", distverif1, distverif2); // getchar();



for(t = 0; t < maxtimesimul; t++){
    // COMPUTE AVG OVER 10y //
    avgtW1 += (avgW1[t])/(maxtimesimul);
    avgtW2 += (avgW2[t])/(maxtimesimul);
    avgtW3 += (avgW3[t])/(maxtimesimul);
    avgtW4 += (avgW4[t])/(maxtimesimul);

    avgtT1 += (avgT1[t])/(maxtimesimul);
    avgtT2 += (avgT2[t])/(maxtimesimul);
    avgtT3 += (avgT3[t])/(maxtimesimul);
    avgtT4 += (avgT4[t])/(maxtimesimul);

    avgtY1 += (avgY1[t])/(maxtimesimul);
    avgtY2 += (avgY2[t])/(maxtimesimul);
    avgtY3 += (avgY3[t])/(maxtimesimul);
    avgtY4 += (avgY4[t])/(maxtimesimul);

    avgtK1 += (avgK1[t])/(maxtimesimul);
    avgtK2 += (avgK2[t])/(maxtimesimul);
    avgtK3 += (avgK3[t])/(maxtimesimul);
    avgtK4 += (avgK4[t])/(maxtimesimul);

    avgtBK1 += (avgBANK1[t])/(maxtimesimul);
    avgtBK2 += (avgBANK2[t])/(maxtimesimul);
    avgtBK3 += (avgBANK3[t])/(maxtimesimul);
    avgtBK4 += (avgBANK4[t])/(maxtimesimul);
    
    avgtMPLDL1 += (avgMPLDL1[t])/(maxtimesimul);
    avgtMPLDL2 += (avgMPLDL2[t])/(maxtimesimul);
    avgtMPLDL3 += (avgMPLDL3[t])/(maxtimesimul);
    avgtMPLDL4 += (avgMPLDL4[t])/(maxtimesimul);

    growthtK1 += (growthK1[t])/(maxtimesimul);
    growthtK2 += (growthK2[t])/(maxtimesimul);
    growthtK3 += (growthK3[t])/(maxtimesimul);
    growthtK4 += (growthK4[t])/(maxtimesimul);
        
    growthtY1 += (growthY1[t])/(maxtimesimul);
    growthtY2 += (growthY2[t])/(maxtimesimul);
    growthtY3 += (growthY3[t])/(maxtimesimul);
    growthtY4 += (growthY4[t])/(maxtimesimul);

    surv1t += (surv1[t])/(maxtimesimul);
    surv2t += (surv2[t])/(maxtimesimul);
    surv3t += (surv3[t])/(maxtimesimul);
    surv4t += (surv4[t])/(maxtimesimul);
}






//////////////////////////////////////////////////////////////////////////////
//////////////////////// WRITE ALL INDICATORS IN FILE ////////////////////////
//////////////////////////////////////////////////////////////////////////////

///////////// PRINT FULL STATISTICS ////////////////
FILE *simulfile;

#if logprint == 2
simulfile=fopen(simul10_2, "w");
setbuf ( simulfile , NULL );


fprintf(simulfile,"surv1\tsurv2\tsurv3\tgY1\tgY2\tgY3\tavgT1\tavgT2\tavgT3\tavgW1\tavgW2\tavgW3\n");
for(t= 0; t<maxtimesimul; t++) {
for(ii = 0; ii < maxgrid; ii++)
{
    for(tt = 0; tt < maxtheta; tt++)
    {
        for(ee = 0; ee < maxfirmtype; ee++)
        {
            fprintf(simulfile, "%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", surv_group1[inxETsimul(ii,tt,ee,t)], surv_group2[inxETsimul(ii,tt,ee,t)], surv_group3[inxETsimul(ii,tt,ee,t)], growthY_group1[inxETsimul(ii,tt,ee,t)], growthY_group2[inxETsimul(ii,tt,ee,t)], growthY_group3[inxETsimul(ii,tt,ee,t)], avgT_g1[inxETsimul(ii,tt,ee,t)], avgT_g2[inxETsimul(ii,tt,ee,t)], avgT_g3[inxETsimul(ii,tt,ee,t)], avgW_g1[inxETsimul(ii,tt,ee,t)], avgW_g2[inxETsimul(ii,tt,ee,t)], avgW_g3[inxETsimul(ii,tt,ee,t)]);
        }
    }
}
}

fclose(simulfile);
#endif
    
    
char buffer_simul_10y[150]; // The filename buffer.
snprintf(buffer_simul_10y, sizeof(buffer_simul_10y), "OUTPUT/statistics_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);

char buffer_survival_10y[150]; // The filename buffer.
snprintf(buffer_survival_10y, sizeof(buffer_survival_10y), "OUTPUT/survival_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);

char buffer_size_10y[150]; // The filename buffer.
snprintf(buffer_size_10y, sizeof(buffer_size_10y), "OUTPUT/size_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);
    
char buffer_skill_10y[150]; // The filename buffer.
snprintf(buffer_skill_10y, sizeof(buffer_skill_10y), "OUTPUT/skill_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);
    
char buffer_wealth_10y[150]; // The filename buffer.
snprintf(buffer_wealth_10y, sizeof(buffer_wealth_10y), "OUTPUT/wealth_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);

char buffer_prod_10y[150]; // The filename buffer.
snprintf(buffer_prod_10y, sizeof(buffer_prod_10y), "OUTPUT/prod_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);

char buffer_bankruptcy_10y[160]; // The filename buffer.
snprintf(buffer_bankruptcy_10y, sizeof(buffer_bankruptcy_10y), "OUTPUT/bankruptcy_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);

char buffer_MPL_10y[150]; // The filename buffer.
snprintf(buffer_MPL_10y, sizeof(buffer_MPL_10y), "OUTPUT/MPL_10y_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_type=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, type_simul, rhostar, pLTpar);
    
    
        
// print survival rates for each group:
simulfile = fopen(buffer_survival_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", surv1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", surv2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", surv3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);
    
simulfile = fopen(buffer_size_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgK1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgK2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgK3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);
    
simulfile = fopen(buffer_skill_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgT1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgT2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgT3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);

simulfile = fopen(buffer_wealth_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgW1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgW2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgW3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);

simulfile = fopen(buffer_prod_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgY1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgY2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgY3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);

simulfile = fopen(buffer_bankruptcy_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgBANK1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgBANK2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgBANK3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);

simulfile = fopen(buffer_MPL_10y, "w");
setbuf ( simulfile , NULL );
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgMPLDL1[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgMPLDL2[t]);} fprintf(simulfile,"\n");
    for(t = 0; t < maxtimesimul; t++){fprintf(simulfile,"%20.15f,", avgMPLDL3[t]);} fprintf(simulfile,"\n");
fclose(simulfile);
    

    
simulfile=fopen(buffer_simul_10y, "w");
setbuf ( simulfile , NULL );

// AVERAGE WEALTH //
fprintf(simulfile,"-------AVERAGE THETA --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", avgtT1, avgtT2, avgtT3, avgtT4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgT1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgT2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgT3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgT4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// AVERAGE WEALTH //
fprintf(simulfile,"-------AVERAGE WEALTH --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", avgtW1, avgtW2, avgtW3, avgtW4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgW1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgW2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgW3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgW4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// AVERAGE Y //
fprintf(simulfile,"-------AVERAGE Y --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", avgtY1, avgtY2, avgtY3, avgtY4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgY1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgY2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgY3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgY4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// AVERAGE CAPITAL //
fprintf(simulfile,"-------AVERAGE K --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", avgtK1, avgtK2, avgtK3, avgtK4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgK1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgK2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgK3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgK4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// AVERAGE MPLDL //
fprintf(simulfile,"-------AVERAGE MPLDL --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", avgtMPLDL1, avgtMPLDL2, avgtMPLDL3, avgtMPLDL4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgMPLDL1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgMPLDL2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgMPLDL3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", avgMPLDL4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// AVERAGE GROWTH K //
fprintf(simulfile,"-------AVERAGE GROWTH K --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", growthtK1, growthtK2, growthtK3, growthtK4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthK1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthK2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthK3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthK4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// AVERAGE GROWTH Y
fprintf(simulfile,"-------AVERAGE GROWTH Y --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", growthtY1, growthtY2, growthtY3, growthtY4);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthY1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthY2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthY3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", growthY4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

// SURVIVAL RATE //
fprintf(simulfile,"-------SURVIVAL RATE --------\n");
fprintf(simulfile,"GROUP1: %20.15f\t, GROUP2: %20.15f\t, GROUP3: %20.15f\t, GROUP4: %20.15f\n", surv1t, surv2t, surv3t, surv4t);
fprintf(simulfile,"GROUP1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f\t", surv1[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP1 T1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv1_byT[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP1 T2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv1_byT[maxtimesimul + t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP1 T3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv1_byT[2*maxtimesimul + t]);
}
fprintf(simulfile,"\n");


fprintf(simulfile,"GROUP2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv2[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2 T1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv2_byT[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2 T2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv2_byT[maxtimesimul + t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP2 T3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv2_byT[2*maxtimesimul + t]);
}
fprintf(simulfile,"\n");


fprintf(simulfile,"GROUP3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv3[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3 T1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv3_byT[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3 T2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv3_byT[maxtimesimul + t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP3 T3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv3_byT[2*maxtimesimul + t]);
}
fprintf(simulfile,"\n");


fprintf(simulfile,"GROUP4: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv4[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4 T1: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv4_byT[t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4 T2: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv4_byT[maxtimesimul + t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"GROUP4 T3: \t");
for(t = 0; t < maxtimesimul; t++) {
    fprintf(simulfile,"%20.15f,", surv4_byT[2*maxtimesimul + t]);
}
fprintf(simulfile,"\n");
fprintf(simulfile,"\n");

fclose(simulfile);










//////////////////////////////////////////////////////////////////////////////
//////////////////////// FREE ALLOCATED POINTERS HERE ////////////////////////
//////////////////////////////////////////////////////////////////////////////

free(zeroE);

free(avgW_g1);
free(avgW_g2);
free(avgW_g3);
free(avgW_g4);
free(dist_group1);
free(dist_group2);
free(dist_group3);
free(dist_group4);
free(avgT_g1);
free(avgT_g2);
free(avgT_g3);
free(avgT_g4);
free(avgY_g1);
free(avgY_g2);
free(avgY_g3);
free(avgY_g4);
free(avgK_g1);
free(avgK_g2);
free(avgK_g3);
free(avgK_g4);
free(avgMPLDL_g1);
free(avgMPLDL_g2);
free(avgMPLDL_g3);
free(avgMPLDL_g4);
free(growthK_group1);
free(growthK_group2);
free(growthK_group3);
free(growthK_group4);
free(growthY_group1);
free(growthY_group2);
free(growthY_group3);
free(growthY_group4);
free(surv_group1);
free(surv_group2);
free(surv_group3);
free(surv_group4);

free(distEgroup1JNEi);
free(distEgroup1JNE);
free(distEgroup1JEi);
free(distEgroup1JE);

free(distEgroup2JNEi);
free(distEgroup2JEi);
free(distEgroup2JNE);
free(distEgroup2JE);

free(distEgroup3JNE);
free(distEgroup3JE);

free(distEgroup4JNE);
free(distEgroup4JE);



}
